!\Description:
!  Compute the eigenvalues of the current upper Hessenberg matrix
!  and the corresponding Ritz estimates given the current residual norm.
!\Arguments
!  RNORM   real scalar.  (INPUT)
!          Residual norm corresponding to the current upper Hessenberg 
!          matrix H.
!  N       Integer.  (INPUT)
!          Size of the matrix H.
!
!  H       complex N by N array.  (INPUT)
!          H contains the current upper Hessenberg matrix.
!  RITZ    complex array of length N.  (OUTPUT)
!          On output, RITZ(1:N) contains the eigenvalues of H.
!
!  BOUNDS  complex array of length N.  (OUTPUT)
!          On output, BOUNDS contains the Ritz estimates associated with
!          the eigenvalues held in RITZ.  This is equal to RNORM 
!          times the last components of the eigenvectors corresponding 
!          to the eigenvalues in RITZ.
!
!  Q       complex N by N array.  (WORKSPACE)
!          Workspace needed to store the eigenvectors of H.
!  IERR    Integer.  (OUTPUT)
!          Error exit flag from clahqr or ctrevc.
!-----------------------------------------------------------------------
subroutine cneigh(rnorm, n, h, ritz, bounds, q, ierr)
    implicit none
    integer ierr, n
    real rnorm
    complex bounds(n), h(n, n), q(n, n), ritz(n)
    ! local variables
    complex, allocatable :: hess(:,:)
    complex vl(1), y(2*n)
    real, external :: scnrm2
    real rwork(n), temp
    logical select(1)
    integer j

    !----------------------------------------------------------%
    ! 1. Compute the eigenvalues, the last components of the   |
    !    corresponding Schur vectors and the full Schur form T |
    !    of the current upper Hessenberg matrix H.             |
    !    zlahqr returns the full Schur form of H               |
    !    in WORKL(1:N**2), and the Schur vectors in q.         |
    !----------------------------------------------------------%
    allocate(hess(n, n))
    call clacpy('All', n, n, h, n, hess, n)
    q(:,:) = 0.0D0
    do j = 1, n
        q(j, j) = 1.0D0
    end do
    call clahqr(.true., .true., n, 1, n, hess, n, ritz, 1, n, q, n, ierr)
    if (ierr /= 0) return

    ! 2. Compute the eigenvectors of the full Schur form T and
    !    apply the Schur vectors to get the corresponding eigenvectors.
    call ctrevc('Right', 'Back', select, n, hess, n, vl, n, q, n, n, n, y, rwork, ierr)
    deallocate(hess)
    if (ierr /= 0) return

    !------------------------------------------------%
    ! Scale the returning eigenvectors so that their |
    ! Euclidean norms are all one. LAPACK subroutine |
    ! ztrevc returns each eigenvector normalized so  |
    ! that the element of largest magnitude has      |
    ! magnitude 1; here the magnitude of a complex   |
    ! number (x,y) is taken to be |x| + |y|.         |
    !------------------------------------------------%
    do j = 1, n
        temp = scnrm2(n, q(1, j), 1)
        q(1:n, j) = q(1:n, j) / temp
    end do
    ! Compute the Ritz estimates
    bounds(:) = q(n, :) * rnorm
end
